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Abstract 

In QCD chiral symmetry is explicitly broken by quark masses, the effect of which 
can be described reliably by chiral perturbation theory. Effects of explicit chiral 
symmetry breaking by the lattice regularisation of the Dirac operator, typically 
parametrised by the residual mass, should be negligible for almost all observables 
if the residual mass of the Dirac operator is much smaller than the quark mass. 
However, maintaining a small residual mass becomes increasingly expensive as 
the quark mass decreases towards the physical value and the continuum limit is 
approached. We investigate the feasibility of using a new approximately chiral 
Dirac operator with a small residual mass as an alternative to overlap and 
domain wall fermions for lattice simulations. Our Dirac operator is constructed 
from a Zolotarev rational approximation for the matrix sign function that is 
optimal for bulk modes of the Hermitian kernel Dirac operator but not for the 
low-lying parts of its spectrum. We test our operator on various 32 3 x 64 lattices, 
comparing the residual mass and the performance of the Hybrid Monte Carlo 
algorithm at a similar lattice spacing and pion mass with a hyperbolic tangent 
operator as used by domain wall fermions. We find that our approximations 
have a significantly smaller residual mass than domain wall fermions at a similar 
computational cost, and still admit topological charge change. 
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1. Introduction 

Lattice QCD is a phenomenologically successful regularisation of QCD which 
can accurately predict experimental observables. It places space-time on a dis- 
crete lattice, and the continuum theory is recovered as three limits are taken, 
the volume to infinity, the lattice spacing to zero, and the quark masses to 
their physical values. The extrapolations in the lattice spacing and quark mass 
are made more difficult, with increased errors, by explicit breaking of chiral 
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symmetry within the lattice discretization of the Dirac operator. The Niclson- 
Ninomiya theorem jl}, which states that it is impossible to represent an odd 
number of quark flavours on the lattice without breaking chiral symmetry is 
avoided for the cheapest commonly used families of lattice actions either by 
simulating additional degenerate fermions (staggered fermions) or by explicitly 
breaking chiral symmetry (Wilson fermions). 

There is, however, an alternative, to use a lattice chiral symmetry [2J, defined 
by a Ginsparg- Wilson relation Q, which has a smooth limit to the continuum 
chiral symmetry. The only known and practical lattice Dirac operators which 
satisfy lattice chiral symmetry are discontinuous, being built from the matrix 
sign function. The simplest of these is the overlap operator, first proposed 
by Neuberger and Narayanan 0, 0]. However, the computational difficulties 
involved in simulating overlap fermions mean that their use in lattice simulations 
remains challenging. A previous solution with a a good approximate chiral 
symmetry, the domain wall fermion, which placed left handed and right handed 
fermions on the four dimensional surfaces at the boundaries of a five dimensional 
lattice [f| . This formalism reduces to a form of the overlap operator for infinite 
fifth dimension, giving exactly chiral fermions. It is also equivalent to the four 
dimensional Kenney-Laub-Neuberger 0, H[ Dirac operator. 

There are, however, many other ways in which an approximate chiral sym- 
metry can be maintained on the lattice. There are no a priori reasons to suspect 
that the historical domain wall fermion should be the best (the 'best' implemen- 
tation can be defined as the one that requires the least computational effort for 
a particular residual mass, the standard measurement of explicit chiral symme- 
try breaking on the lattice). In this work, we propose and investigate a four 
dimensional approximation to the overlap operator that gives a considerably 
better residual mass than domain wall fermions while avoiding the complexities 
of exact overlap fermions. 

1.1. On- Shell Chiral Symmetry 

For a lattice Dirac operator lp, an on-shell chiral transformation can be 
written as |2j 

^ ^ e «*75(l-a-0)^ and ^^.^(1-^)75. 

the differential condition that the Dirac operator itself is invariant under this 
transformation is a Ginsparg- Wilson relation 3] 

{ 7 5,0} = 2a0 7 50, 
which may be written in the equivalent forms 

Relp^a^i) <=> Re[(l - a^a$] =0 <^ Re0 1 = a. 

We may define the quantity 75 by alp = 5(1 + 7575), and if we wish lp to satisfy 
lp = 7 50 7 5 then 75 =<%■ This Ginsparg- Wilson relation implies that 7 | = 1, 
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so 75 is both hermitian and unitary, so that its spectrum can contain only the 
eigenvalues ±1 and hence 75 = sgn(7 5 ). 

We must also require that ij) is a Dirac operator in the naive continuum limit, 
Ip = Zqw{$ + A ) + 0(a 2 ), where there are no corrections of 0(a) because the 
only available such operator is [0 + Jft ,0 +4-} = <r ■ F = a a ^F a ^ which is not 
chirally invariant. If Ipw is any 75-hermitian lattice Dirac operator satisfying 
Ipw — Zw{0 +4-) + 0(a) then we may choose 2aMlfi — lp w + O(a) where 
2a M = Zw/Zqw is a suitable finite wave- function renormalization, and we 
obtain the (massless) Neuberger operator [H, Q 

fl o = i(l+ 75 sgnff) (1) 

where H = j$lfiw ~ M- In QCD a quark has a (small) mass /i, and this can be 
incorporated by shifting the spectrum to lie in the interval [a[i, 1] by defining 
the massive Neuberger operator to be 

a&n = |[(1 + om) + (1 - aM)75Sgn(F)]. (2) 

Within this framework there are numerous choices that have to be made [lol | : 
whether to use four or five dimensional pseudofermions; the choice of the ker- 
nel operator H (including the mass M); the choice of rational approximation 
for the matrix sign function; which form of 5D matrix with the desired Schur 
complement (for example, Euclidean-Caley or continued fraction); and for four 
dimensional pseudofermions whether to use a nested 4D or 5D inverter. All 
of these choices are independent and physically equivalent. We use the name 
"domain wall operator" as any approach which uses five dimensional pseud- 
ofermions, and "overlap operator" as any approach in four dimensions with an 
exact matrix sign function, limited only by the floating point precision of the 
computer. An "approximate overlap operator" is a four dimensional approach 
which has a small explicit breaking of chiral symmetry from the use of an ap- 
proximate matrix sign function. In this article, we are only studying the effect of 
changing the rational approximation used for the matrix sign function. Our re- 
sults should be independent of the choice of kernel and whether four dimensional 
or five dimensional pseudofermions are used. 

1.2. Ginsparg-Wilson Defect and Residual Mass 

Suppose we have some approximation s(H) 1=3 sgn(if) to the matrix sign 
function, so that equation ([2]) is replaced by 

D li = ±[(l + an) + (l-a l i)e(H)] ) (3) 

then we may define the defect A as the amount by which the corresponding ap- 
proximate Neuberger operator, J]) , fails to satisfy the Ginsparg-Wilson relation 
for /i = 0, 

A ee |{75,-0o} - a$o750o = 75 Re[(l - a$^$ Q ], 
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which gives 

4 75 aA = 1 -e(H) 2 . (4) 
From equation © we have $ ^ = $0 + — alfio), hence 

Re[(l - alpo)^,.} = Rc[(l- 000)^0] +M(l-a0o) t (l-a0o) 
= 7 5 A + /i(l-a^ o ) t (l-a0o)- 



Multiplying this by l]) t 1 on the left and J]) „ 1 on the right we obtain 

±(S M + St) = ^t-^sA^- 1 + M Sj5„; (5) 

where Sp = (1 — alpo)^^ 1 (which satisfies S^ 1 = Sq + fi). The first term on 
the right side of equation ([5]) is a measure of the chiral symmetry breaking due 
to the approximation to the sign function whereas the second is that due to the 
explicit quark mass afi. We wish to introduce some norm ||A|| on the defect 
to quantify the magnitude of the errors due to our approximation to the sign 
function. A useful estimate for this norm is the residual mass 

, tr^t-^A^-i) 

TTL — — 

tr(StS M ) 

whence tr(S M + St)/2 tr(S M St) = m ros + fi. 

There is no a priori reason when comparing the extent of chiral symmetry 
breaking for different Dirac operators why we should project the defect onto 
a scalar using this trace, indeed in our numerical studies we have used the 
computationally cheaper approach of projecting onto momentum zero states, 
giving 
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We do not expect that the quantity m ros will be less suitable than m' icB to com- 
pare the effects of chiral symmetry breaking between different Dirac operators. 

Since our approximation for the Neuberger operator does not exactly satisfy 
chiral symmetry it will have 0(a) corrections near the continuum limit, and as 
these must either come from the explicit mass term or the defect we see that 

$V = Z GW (0 + 4 ) + M + m xes + constant x aa ■ F + 0(a 2 ). (6) 



1.3. Effects of Residual Mass 

The spectrum of the exact unitary matrix 7575 = 75 sgn(iJ) lies on the 
unit circle in the complex plane 1 1 TsO's 1 1 = 1. The approximate matrix sign 
function has spectral norm = sup|^,| =1 ip^£(H)ip, so the spectrum of 

||75£(-ff)|| < INI \\£(H)\\ = \\e(H)\\ lies within a disc of radius \\e(H)\\. For 
an operator with good chiral symmetry most eigenvalues will be close to the 
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unit circle; however small eigenvalues of H, where the approximation is less 
good, may exhibit large discrepancies from the spectrum of the exact overlap 
operator. For example, if H has an exactly zero eigenvalue then for all symmetric 
approximations to the matrix sign function the corresponding eigenvalue of Ij) ^ 
will lie in the centre of the circle, i.e. at |(1 + a[i). 

If the residual mass is sufficiently smaller than the target physical quark 
mass, m roa <C /J. then as long as the lattice spacing is small enough that the 0(a 2 ) 
effects in equation ([6]) are negligible all physical effects of m tes can be removed 
by adjusting \i h-> fx — m roa so that the quark mass stays fixed. However, if the 
residual mass is larger than the physical quark mass then this is not possible in 
general. While one could insert \i < into the lattice Dirac operator this would 
be likely to introduce zero or negative eigenvalues, leading to the inexact overlap 
Dirac operator becoming singular on some (now exceptional) configurations. 

It is therefore desirable to have a small residual mass. The difficulty with 
this is that equation (j4]) tells us that ||A|| = ||1 — e(H)\\/4a, so we need to reduce 
the error || 1 — e(H ) || oc a in our approximation to keep the physical m roa fixed as 
the continuum limit is approached. This can be done by increasing the order of 
the rational approximation which may, depending on the approximation used, 
significantly increase the cost of the simulation. 

At larger lattice spacings it is easier to maintain a small m roa , but there will 
be larger 0(a 2 ) and higher lattice artefacts in equation ([6]); since these come 
from the chiral symmetry breaking in A and not just from the quark mass /i 
they are less easy to model using, for example, chiral perturbation theory. 

1.4- Choice of Rational Approximation 

We have shown that we need a lattice Dirac operator with a good approx- 
imation to (on-shell) chiral symmetry, but we also wish to avoid the cost of 
maintaining exact chiral symmetry to machine precision. Our goal is thus to 
find a family of approximations to the matrix sign function that provides a good 
balance between residual mass and the cost of Hybrid Monte Carlo (HMC) com- 
putations. We shall demonstrate that domain wall fermions, where the accuracy 
of whose hyperbolic tangent approximation falls too slowly with increasing ra- 
tional degree are far from optimal, while overlap fermions, which are far more 
chiral than necessary for massive quarks at a cost in time and complexity of the 
HMC algorithm, are more expensive than is required to obtain a good enough 
chiral symmetry. 

In this work, we introduce and test the Zolotarev lattice Dirac operator, 
which uses the optimal rational approximation to the matrix sign function, but 
not over the entire spectrum of the kernel operator ifi w ■ This is guaranteed to 
provide the best approximation to the matrix sign function for a given order of 
rational approximation within a certain tunable eigenvalue range, and therefore 
might be expected to give the smallest residual mass for a given amount of 
computational effort. However, there are a number of questions which need to 
be addressed in this comparison. 
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1. While the Zolotarev approximation is guaranteed to give the best L m 
approximation over part of the spectrum, some eigenvalues lie outside 
this range. What contribution do they make to the residual mass? 

2. Does the reduced residual mass come at the cost of less stable molecular 
dynamics? If a small time step were required to resolve larger fluctuations 
in the fermionic force this could compensate for any gain in the residual 
mass. 

3. Is the Zolotarev Dirac operator ergodic? In particular, are all topological 
sectors sampled, and is the autocorrelation between topological sectors 
and different portions of the same topological sectors short enough? 

4. Is the Zolotarev Dirac operator local? As it does not approximate the sign 
function well for small eigenvalues of the kernel operator it is not obvious 
that we can rely on previous proofs and numerical results, although it 
seems unlikely that it could be less local than the hyperbolic tangent 
Dirac operator used in domain wall computations. 

In [J3]we describe our Zolotarev approximation and compare it to the Kenney- 
Laub-Neuberger @,[B[ hyperbolic tangent approximation, which is a four dimen- 
sional representation of the five dimensional domain wall Dirac operator; in £[3] 
we describe how we tested the various methods; in 21 we present the numerical 
results from these tests; and in Sj5]we present our conclusions. 

2. The Zolotarev Dirac operator 

2.1. The Zolotarev approximation 

The Zolotarev Dirac operator depends on the degree N of the rational ap- 
proximation and a parameter £ where £,\\H\\ < |A| < \\H\\ is the interval of the 
kernel operator's spectrum on which the approximation to the matrix sign func- 
tion is optimal. Here \\H\\ = sup (u, Hu) is the spectral norm of the hermitian 

IM|2 = 1 _ 

kernel operator i.e., its largest eigenvalue. The Zolotarev Dirac operator is 



where 



(1 + afx) + (1 - a^)75 e 2 



H 

W\ 



e z (x) =x 



L7V/2J 



x 2 - af 

2 — 1 1 



is the Zolotarev approximation to the sign function, which is optimal in the L m 
norm over £ < |x| < 1. The coefficients Wj and Cj can be computed in terms of 
Jacobi elliptic functional depending upon £ and N [Tl, 12, 13, 14 1. Were the 



1 In particular Uj = £sn^2iK'(j — j)/N t ^\, where K'(fc) = K (Vl — k' 2 ^ is a complete 
elliptic integral. 
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Figure 1: The solid lines show Zolotarcv rational approximations Ez{x) for various degrees TV 
with £ = 10~~ 2 as a function of log 10 x. For comparison the dashed lines show the Kenney- 
Laub- Neuberger hyperbolic tangent approximations ekln (x) of the same degrees. 



Zolotarev rational approximation expressed in a five dimensional formalism, 
similar to that used for domain wall fermions, the size of the fifth dimension 
would be L s = N. 

There are two families of Zolotarev approximations, those that vanish at 
the origin and those that are singular at the origin (up to a small re-scaling 
the latter are just the reciprocals of the former). We shall only use the non- 
singular kind, which vary smoothly and monotonically from —1 just below — £ 
to +1 just above £. Figure [T] shows the Zolotarev approximation e z (x) for 
various degrees N and for a typical value £ = 0.01. The maximum error of 

e z over the interval on which it is optimal, A = max |e z (x) — 1|, is shown in 

£<W<i 

Figure [2] Not only can we see from Figure [2] that the error in the approximation 
A falls exponentially with the degree N over the interval where the Zolotarev 
approximation is optimal, but also from Figure [1] we see that the approximation 
improves over the interval \x\ < £. 



2. 2. The KLN Dirac operator 

The Kenney-Laub-Neuberger Dirac operator operator (KLN), which is im- 
plicitly used in the domain wall approach 1^, 15, 16 1 , is 



(1 + a/x) + (1 - a/x)75 



H 

W\ 
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Figure 2: The maximum error of the Zolotarev approximation Z(x) over its approximation 

interval, A = max |ez(a;) — 1| for various degrees N are shown on a log-log plot. 
f<kl<l 



where the Kenney-Laub-Neuberger 0, [B| hyperbolic tangent approximation to 
the sign function is 

£ KLN (z) = tanh (N tanh- 1 x) = ^i±|l^— 1L-|L , 

This may be written for even N as a partial fraction expansion (lpj 



2x i-^ 



1 



tan 



iV 



V S*»+(tan^y 



Figures Q] shows how this approximation compares to the Zolotarev approx- 
imation. Unlike the Zolotarev approximation it does not require a minimum 
eigenvalue as an input; like the Zolotarev approximation the accuracy depends 
on the order of the rational approximation. 

2. 3. Fermionic forces in the Hybrid Monte Carlo algorithm 

One of the major possible difficulties with this method is that a very small 
integration step size St may be required to prevent the molecular dynamics 
(MD) trajectory in the Hybrid Monte Carlo (HMC) il?} algorithm from be- 
coming unstable. For overlap fermions these instabilities can be avoided (albeit 
with a little additional complexity in the force calculation [18J) because the 
small eigenvalues are not treated by an approximation but deflated, and the 
simulation is run in the chiral sector without zero modes [HH^. Here, because 
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all the eigenvalues are treated by the rational approximation, the force acting 
on the gauge fields from their interaction with fcrmion fields might become large 
for three reasons: 

1. the fermion mass is small and the gauge-pseudofermion coupling involves 
the inverse of the Dirac operator which will, in general, have approximate 
zero modes; 

2. the derivative of the chiral Dirac operator ij) ^ becomes large if the fermion 
field is close to a zero mode of the kernel operator H (for an exact overlap 
operator, it is a Dirac S- function); and 

3. the estimate of the fermionic force obtained from pseudofermion fields is 
noisy 



21|,|22| 



The use of multiple pseudofermion fields [23j, |24| can reduce the effects of 



and, for overlap fermions, a factorisation of the determinant can can completely 



remove the effect of the pseudofermion noise [22J. In overlap simulations, various 
transmission-reflection methods 2|| 2^, [2(| have been introduced to resolve the 



Dirac ^-function from the effects of point [2] above. However, the adaptations to 
the usual HMC algorithm needed for exact overlap fermions are expensive, so 
for the case of interest here — where chiral symmetry is broken explicitly by a 
small fermion mass afi — we advocate choosing a Zolotarev approximation that 
is not optimal for the smallest eigenvalues of H, but is a compromise between 
having a small residual mass m ros <C /i and having a small fermionic force. It is 
clear from Figure [1] that we may expect the Zolotarev Dirac operator to have a 
much smaller residual mass than the KLN Dirac operator used in the domain 
wall method for any degree N and reasonable values of £. 

The fermionic force for continuous time MD evolution, to which the discrete 
integrators used in HMC are a good approximation for reasonable acceptance 
rates, contains a term proportional to the derivative of the approximation to the 
sign function used in [18]. The derivatives of the Zolotarev and KLN approx- 
imations are shown in Figure [3J and several interesting features arc immediately 
obvious. For both Zolotarev and KLN approximations the largest derivative oc- 
curs at x — 0, this is trivial for the KLN approximation where the derivative 
at the origin is £' KljN (0) — N, but it is less obvious for the Zolotarev approxi- 
mation. We first note that the derivative e' z (x) over the interval £ < |x| < 1 
on which the approximation is optimal is always small, and falls exponentially 
with the degree N. Next we see that the derivative increases monotonically as 
x approaches zero where it attains its maximum value. In Figure [4] we show the 
dependence of e' z (0) on £ and the degree N; empirically the data are well fitted 
by 4(0) « (0.49 + 0.025iV)/£ - 87 . 



3. Numerical Implementation 

3.1. Markov Chains 

Our goal in this work is to investigate whether the Zolotarev Dirac operator 
gives a smaller residual mass for the same cost as the domain wall Dirac opera- 
tor. To do this, we have generated several small ensembles on 32 3 x 64 lattices 
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Figure 3: The derivative of the Zolotarev rational approximation \de z {x)/dx\ with § = 10 — 2 
and of the Kenney-Laub- Neuberger hyperbolic tangent approximation \deKLN(x)/dx\ for var- 
ious orders N. 




Figure 4: The derivative of the Zolotarev rational approximation \dez/dx\ at the origin x = 
for various orders N . 
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starting from equilibrated domain wall configurations. We aimed to match the 
pion masses and lattice spacings for these ensembles, although in practice the 
Zolotarev ensembles were at a slightly lighter pion mass and slightly smaller 
lattice spacing. We compared the KLN Dirac operator, which is equivalent to a 
domain wall fermion with a Borici kernel, with three different Zolotarev Dirac 
operators with different values of £. We used a Borici- Wilson kernel with one 
step of over-improved [27| stout smearing 
proved Liischer-Weisz gauge action 29|, 30, 



at k = 0.19 and a tadpole im- 
Q . We ran enough trajectories 
in each case to ensure that the plaquette was thermalised with the new action 
before taking measurements of the pion mass, lattice spacing, and residual mass. 
There is a small systematic uncertainty in these measurements as our ensembles 
were not fully thermalized with respect to the measured observables, however 
we have noticed no change in our measurements along the Markov chains, so 
any systematic error due to incomplete thermalization is probably smaller than 
our statistical errors. We have made measurements on 5-10 configurations for 
each run. Our goal in making these measurements is not to extract physics, but 
to obtain an approximate idea of the physical parameters. 



3.2. MD Integrators 

Construction of the HMC MD integrator was straightforward. On the larger 
lattices, we used three additional pseudofcrmions, following the method of Hasen- 
busch 23] and multiple time-scale integration [33j with the gauge field time steps 
being eight times smaller than that for the the heaviest pseudofermion, which in 
turn had a time step eight times smaller than that of the lightest pseudofermion. 



3.3. Linear Equation Solvers 

For any choice of rational approximation and kernel operator there are sev- 
eral different approaches to the problem of inverting the chiral Dirac operator 
JP^ 0. 

1. Introduce a five dimensional matrix that has Ip ^ as its Schur complement. 

(a) Use this to find the inverse of Ip ^ applied to a four dimensional 
pseudofermion field in order to compute the fermionic force for four 
dimensional MD. 

(b) Use the five dimensional matrix as part of a five dimensional MD 
scheme with five dimensional pseudofermions (this corresponds to 
the domain wall formalism). 

2. Apply the inverse of Ip ^ to a four dimensional pseudofermion field by use 
of a nested solver. 

In this paper we use a nested four dimensional solver, using a par tial fraction 
representation with relaxation [34| . GMRESR preconditioning |35j |. and defla- 
tion of about 70 kernel eigenvalues. 

It is not in the remit of this paper to compare the performance of four- 
and five-dimensional solvers; our goal is to compare the effect of different ratio- 
nal approximations. We expect that the relative performance of such different 
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Volume 




io 4 e 


N 


a (fm) 




k 


Acc. % 


32 3 x 64 


8.65 


KLN 


16 


0.101(3) 


0.203(3) 




90 


32 3 x 64 


8.7 


1.93 


16 


0.094(4) 


0.182(5) 


35 


100 


32 3 x 64 


8.7 


0.643 


16 


0.096(3) 


0.182(1) 


18 


97 


32 3 x 64 


8.7 


0.193 


16 


0.095(3) 


0.183(2) 


9 


92 



Table 1: Parameter values for ensembles, p is the Liischer- Weisz gauge coupling, a and m w 
are the lattice spacing and pion mass respectively. We also show N, the degree of the rational 
approximation, and the HMC acceptance rate. For the Zolotarev Dirac operator we list £ and 
k, the average number of eigenvalues of the kernel operator H that lie outside the optimal 
interval of the Zolotarev approximation. The first line corresponds to our KLN approximation 
ensembles. 

approximations to the sign function should be similar with the use of five di- 
mensional solvers such as used in the domain wall algorithm. 

The parameters of the runs are given in Table [T] Pion masses were measured 
using the pseudo-scalar correlator, and the lattice spacing using ro [36l l37j |. 
Because we only have a handful of configurations for each ensemble, it was 
difficult to get a reliable estimate of the statistical errors, particularly for the 
lattice spacing, and a larger study is required to get more accurate values. 
However, the physical pion masses, which are in the range 370 — 400 MeV, and 
the lattice spacings agree across our ensembles to within the errors. 

4. Results 

Jf..l. HMC Acceptance Rate 

The change in the value of the Hamiltonian AE for our HMC runs is shown 
in Figure [5j There are no large spikes, nor any large differences between the 
various different Dirac operators. The acceptance rate shown in Tables [T] are 
similar (albeit high) for each of the HMC runs. The ensembles used the same 
Hasenbusch masses, MD time steps and, for our performance tests, the same 
number of deflated eigenvalues. There is no indication of any MD instabilities. 

4-2. Performance 

The total time taken for the inversions required for each HMC run for the 
32 3 x 64 lattices is shown in Table [2] The number of projected Wilson eigen- 
values varied from trajectory to trajectory, with an average around 70. All the 
measurements were performed with identical time-steps, trajectory lengths, and 
an equal number of pseudofermion fields. While generating the configurations, 
the order of all the Zolotarev Dirac and KLN operators was held fixed at 16. 

While a larger study is needed for definitive results, it is clear that the 
computational cost required for the Zolotarev and KLN runs does not differ by 
any substantial amount. The observed, small, difference is caused by a slower 
convergence for the inversion of the Zolotarev Dirac operators compared to the 
hyperbolic tangent operator, which is partially explained by the slightly heavier 
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0.40 



■ Zolotarev 0000193 

□ Zolotarev 0.0000643 

□ Zolotarev 0.0001 93 

□ KLN 




-0.47 -O.40 -0.33 -0.27 -0.20 -0.13 -0.07 0.00 0.07 0.13 0.20 0.27 0.33 0.40 0.47 
AE 



Figure 5: The distribution of the HMC energy difference for the KLN and three Zolotarev 
Dirac operators for the 32 3 X 64 lattice. 



Volume 


10 4 £ 


Inverter 


Eigenvalues 


32 3 x 64 


KLN 


5,435 


1,568 


32 3 x 64 


1.93 


6,291 


1,627 


32 3 x 64 


0.643 


6,650 


2,657 


32 3 x 64 


0.193 


6,194 


1,641 



Table 2: Time (in seconds) spent in the two dominant parts of the HMC code in each trajec- 
tory. 
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• Zolotarev 0.0000193 -•- Zolotarev 0.0000643 V Zolotarev 0001 93 -*-KLN 
0.25-1 




Figure 6: The distribution of the smallest eight eigenvalues of the kernel operator H during 
the molecular dynamics for the KLN and Zolotarev Dirac operators. 

quark mass for the KLN. From these studies we conclude that the difference in 
cost between running Zolotarev and KLN fermions at equal degree and equal 
pion mass is at most about 20%. We must stress that these inversions were 
not performed at equal residual mass; to run at the same residual mass as the 
Zolotarev fermions, the cost for the KLN or domain wall fermions would be far 
greater (see §4.5|) . 

4-3. Rate of topological charge changes 

In Figure [H we plot the distribution of the ten smallest eigenvalues of the 
kernel operator H as it evolved during the molecular dynamics. From this lim- 
ited data we cannot reliably estimate the rate of topological charge change. 
The large fluctuations in the data are due to our small sample size; nonethe- 
less the data is clear enough for a qualitative picture to emerge. Our data is 
consistent with our expectation that for the Zolotarev with largest £ there is 
no suppression of the small kernel eigenvalues compared with KLN; although 
there are indications that the eigenvalues of the Zolotarev with the smallest £ 
might be suppressed. No difference can be seen from our data between the rate 
of topological tunnelling for the Zolotarev at large £ and KLN Dirac operators. 

The suppression of the small eigenvalues also aids the stability of the HMC 
algorithm and the locality of the Zolotarev operator. 

4- 4- Locality of Dirac operators 

We checked the locality of the Dirac operator by measuring the exponential 
decay from a single source. All the Dirac operators are exponentially local, and 
there is no noticeable difference in locality between the overlap operator and 
the Dirac operators studied here. 



14 



I Zolotarev 0.00001 93 ♦ Zolotarev 0.0000643 V Zolotarev 0.0001 93 A KLN 




Figure 7: The residual mass m ro3 plotted as a function of the degree N of the KLN and three 
Zolotarev approximations, as measured on the 32 3 X 64 lattices. 



4-5. Chiral symmetry breaking 

The residual mass m res , as denned in H1.2I is shown for both the KLN and 
all three Zolotarev Dirac operators as a function of the degree of the rational 
approximation in Figure [7] In this plot, m rcs has been averaged over all con- 
figurations. It is immediately seen that both Zolotarev Dirac operators exhibit 
significantly smaller residual masses than the KLN operator for all but the low- 
est degrees. At the order used in our tests, N = 16, equivalent to Ls = 16 
for domain wall fermions, the improvement is four orders of magnitude. We 
also observe a striking difference in the way in which the chiral symmetry is 
broken for the KLN and the Zolotarev Dirac operators. For the KLN Dirac 
operator the residual mass decreases steadily. For the Zolotarev Dirac operator 
it decreases rapidly up to around degree 16, and, for the larger Zolotarev ranges 
reaches a plateau, where the fluctuations are dominated by statistical noise, and 
for the smallest range continues to decrease but at a slower rate as the order 
is increased. Moreover, while m ros is roughly constant between configurations 
for the KLN Dirac operator, there is a large fluctuation between configurations 
for the Zolotarev Dirac operator. There is therefore little point in running a 
Zolotarev approximation, at least for the parameter values considered here, with 
degrees N > 16. 

The pattern of chiral symmetry breaking for the Zolotarev Dirac operator 
can be easily explained: there are two contributions to the violation of the 
Ginsparg- Wilson relation for the Zolotarev Dirac operator, from the imper- 
fection of the approximation to the matrix sign function within the interval on 
which the Zolotarev approximation is optimal (the "bulk"), and the contribution 
from the small eigenvalues below £ \\H\\. The first contribution decreases rapidly 
with TV, while the second contribution decreases considerably more slowly, if at 
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Figure 8: The ratio of the contributions to m roB from bulk modes to that from low modes 
as a function of the degree of the rational function approximation N for the KLN and three 
Zolotarev Dirac operators, measured on the 32 3 X 64 lattices. 



all. For small N m ICS is dominated by the bulk, while for large N it is dominated 
by the low modes. 

These two sources of explicit chiral symmetry violation are illustrated in 
Figure El in which we have separated the contributions to m ros from the bulk 
modes from the modes below £\\H\\ and plotted their ratio, defined as 



s = 



V , H>1 75 A(£ i hfc>«'i|)#; 



(7) 



for for the lowest orthonormal eigenvectors, ipi, of the kernel operator H. It can 
be seen that for the Zolotarev approximation for order N > 16 the chiral sym- 
metry breaking is almost entirely caused by the lowest eigenvalues of the kernel 
operator, while for lower N the bulk eigenvalues have a larger contribution. It 
is thus clear that the plateau in m ros for the £\\H\\ = 0.0000643 ensemble is 
caused by the low modes, whereas the S,\\H\\ = 0.0000193 ensemble has no such 
small eigenvalues, and correspondingly the residual mass gradually improves as 
the order of the rational approximation improves and the approximation for the 
matrix sign function gets better for the bulk modes. 



5. Conclusions 

Our main conclusion is that the Zolotarev Dirac operator seems to provide 
a significant improvement over both domain wall and overlap fermions for com- 
putations with light fermions on fine lattices. 
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On coarse lattices contributions from the small eigenvalues subspace of the 
kernel operator dominate the residual mass. Some of these eigenvectors are 
topological, whereas others are lattice artefacts ( "defects" ) . Various mechanisms 
have been suggested to suppress such defects without affecting the underlying 
continuum physics by increasing autocorrelations for topology change. To the 
extent that such methods are effective the contributions to the residual mass 
from the bulk will still be significantly reduced by the Zolotarev Dirac operator. 

5.1. Comparison with Domain Wall fermions 

The principal advantage that the Zolotarev Dirac operator has over domain 
wall fermions is that it gives a significantly smaller residual mass for the same 
computational effort. In particular, this allows for simulations at smaller lattice 
spacing than are currently possible with domain wall fermions. The second 
advantage is that the Zolotarev Dirac operator can be tuned to balance the 
performance of the HMC, the tunnelling rate and the residual mass. 

There are still questions which need to be addressed in future work. In 
particular, we have not addressed the question of whether using five dimensional 
(as used by Chiu and collaborators |16| ) or four dimensional pseudofermions as 
in this work is superior. Should the five dimensional inversion prove superior to 
the nested four dimensional inversion (the comparison in [38} used a considerably 
sub-optimal nested 4D algorithm, so this question remains open), then this can 
easily incorporated into the four dimensional algorithm, so it seems unlikely that 
there is much difference between the two. However, until a direct comparison 
is made in a future work, no definite statement can be made in favour of either 
formulation. 



5.2. Comparison with Overlap fermions 

The Zolotarev Dirac operator has the advantage over overlap fermions that 
it is faster. An exact overlap calculation requires an accurate resolution of the 
matrix sign function during the molecular dynamics, which in practice requires 
the transmission/reflection algorithm [25, 26, 2(j. To allow frequent topological 
charge changes, the overlap HMC algorithm has to be further refined and care- 
fully tuned [22|, leading to approximately a doubling of the cost per trajectory. 
Furthermore, while these refinements allow topological charge changes every few 
trajectories, there may still be longer autocorrelations than with the Zolotarev 
Dirac operator. This cost can be removed by adding unphysical terms to the 
action [39|, thus forbidding both topological tunnelling and kernel eigenvalues 
close to zero. The effect of this with regards to possible artefacts and ergodicity 
is, however, unclear. 
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